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INTRODUCTION 


The  search  tor  "optimal”  shell  finite  elements  has  been  underway  for  nearly  two  decades. 
In  recent  years,  it  has  further  accelerated  in  light  of  significant  progress  in  the  technology  of 
shear-deformable  C°  bending  elements  (e.g..  References  1-19).  Although  the  main  obstacles 
for  these  developments,  known  as  shear  and  membrane  locking  phenomena,  have  been 
addressed  extensively  and  several  remedial  schemes  have  been  proposed,  a  viable  three-node 
doubly  curved  shear-deformable  element,  which  is  the  most  desirable  element  for  general  shell 
analysis,  has  not  yet  been  developed.  The  purpose  of  this  effort  is  to  derive  such  an  element 


We  base  our  finite  element  derivation  upon  Reissner-Mindlin  plate  theory  which  will  con¬ 
stitute  the  bending  part  of  the  element.  To  account  for  the  membrane  deformations  and  the 
membrane-bending  coupling  associated  with  the  shell-element  curvatures,  we  shall  resort  to 
Marguerre's  shallow  shell  equations.  Shallow  shell  elements  of  this  type  specialized  to  the 
axisymmetric  response  proved  effective  in  discretizing  shallow  as  well  as  deep  shell  struc¬ 
tures.1"  The  major  advantage  of  this  analytic  approach  over  general  shell  formulations  (e.g., 
References  5  and  18)  is  its  inherent  simplicity.  Herein,  the  displacements  and  stress  resul¬ 
tants  are  attributed  to  the  element  reference  plane.  Consequently,  integrations  are  carried 
out  across  the  reference  plane  rather  than  the  curved  surface  as  in  the  general  shell  elements. 

According  to  Reissner-Mindlin  theory,20'22  the  strain-displacement  relations  can  be 
expressed  as: 


*  -  <*„• 

if  ,  k  }T  =  lx(6) 

yy  xy  1  ' 
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where  k  and  -y  are,  respectively,  the  curvature  and  transverse  shear  strain  vectors,  8  is  the 
bending  rotation  vector 


0T  ”  {0y’0x} 


(3) 


with  0X  and  0y  denoting  the  bending  rotations  about  the  x  and  y  axes,  respectively,  w  is  the 
transverse  displacement  (refer  to  Figure  1),  and  the  superscript  T  denotes  transpose;  L(  and 
L2  are  the  linear  differential  operators,  and  I  is  an  identity  matrix: 
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The  Marguerrc  membrane  strain-displacement  relations  for  a  thin  shallow  shell  have  the 


form 


23 


e  =  Lx(u)  +  LiU)L2(w) 


(5) 
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where  u  and  v  are  the  membrane  displacements  in  the  x  and  y  coordinate  directions,  respec¬ 
tively;  and  £  =  §  (x,  y)  is  the  initial  height  of  the  shallow  shell. 


X 


Figure  1 .  Shallow  shell  notation. 

One  important  aspect,  which  in  previous  attempts  to  merge  Reissner-Mindlin  and 
Marguerre  theories  has  not  been  addressed,9'12  is  the  conceptual  difference  in  the  transverse 
displacement  variables  appearing  in  Equations  2  and  5.  In  Equation  2,  w  is  a  weighted  aver¬ 
age  transverse  displacement  across  the  thickness,  whereas  in  Equation  5,  w  represents  the  mid¬ 
surface  transverse  displacement.  The  former  variable  comes  into  play  due  to  the  inclusion  of 
shear  deformation  in  Reissner-Mindlin  theory;  the  latter  one  is  a  consequence  of  the  Kirch- 
hoff  thin-regime  assumption,  which  neglects  shear  deformation.  Utilizing  Equation  2,  the 
Kirchhoff  thinness  constraint  reads; 

1*2 (w)  *  -  10  .  (7) 

Replacing  Equation  7  into  Equation  5  yields  the  Marguerre  membrane  strains  consistent  with 
the  Reissner-Mindlin  strains: 

e  —  I*i(u)  ■  (S) 

The  stress  resultants,  which  are  attributed  to  the  reference  plane  of  the  shell,  are  related 
to  the  strains  through  the  constitutive  law: 

N  ■ {N**-  V  vT  *  ^ 
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where  A,  D,  and  G  are,  respectively,"  the  membrane,  bending,  and  transverse  shear  constitutive 
matrices.  For  an  isotropic  shell  of  constant  thickness,  t,  the  constitutive  matrices  arc  given  bv 
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where  E  is  Young’s  modulus,  G  is  the  shear  modulus,  and  v  is  Poisson’s  ratio;  k"  =  5/6  is 
Reissner’s  shear  correction  factor.20 


The  principle  of  virtual  work  can  then  be  employed  to  derive  the  finite  element  stiffness 
equilibrium  equations: 


JJ  (N^5e  +  M^6ic  +  Q^Sy  ~  q6w)  dA  *  0 


(12) 


where  q  is  the  distributed  transverse  loadmg,  A  is  the  reference  plane  area,  and  <5  denotes 
the  variational  operator. 


FINITE  ELEMENT  ISSUES 

The  development  of  effective  curved  shear-deformable  shell  elements  is  severely  hampered 
by  the  “locking  phenomena”  (extreme  stiffening),  reflecting  the  inability  of  the  shell  to  bend 
without  stretching  (“membrane  locking”)  and  transverse  shearing  (“shear  locking”).  The  two 
phenomena  are  directly  linked  to  the  penalized  strain  energy  which,  in  its  nondimensional 
form,  can  be  expressed  as: 

U(ic,y,e)  *  U,  (k)  +  a  U  (y)  +  a  U  (e)  (13) 

D  s  s  m  ID 

in  which  Ub(ic),  Us(y),  and  Um(£)  denote  the  nondimensional  bending,  transverse  shear,  and 
membrane  energy  integrals;  and  as  and  am  are  the  nondimensional  shear  and  membrane  pen¬ 
alty  parameters,  respectively.  Note  that  as  =  0(A2/t2)  and  am  =  0[(x-;A2)2/t"l,  where  A  and  x\- 
are,  respectively,  some  characteristic  span  and  curvature  of  the  shallow  shell.'1'1"  As  the  shell 
thickness,  t,  diminishes  to  zero,  both  as  and  am  approach  infinity,  thereby  enforcing  the  vanish¬ 
ing  shear  and  membrane  strains: 

L2(w)  -+  -10  (Kirchhoff  constraints)  (14a) 

Li(u)  -*•  Lx(C)e  (Membrane  inextensibility  constraints)  .  (14b) 

The  particular  appeal  of  this  theory  is  that  the  variational  statement.  Equation  12, 
requires  a  class  of  C°  continuous  approximations  for  the  w,  u,  and  0  fields  (since  their  high¬ 
est  spatial  derivative  in  Equation  12  is  of  order  one)  and,  therefore,  simple  shape  functions 
can  be  used.  On  the  other  hand,  constraints.  Equation  14,  when  imposed  at  the  element 

level,  pose  severe  limitations  on  the  kinematic  freedom  attainable  by  each  element,  often 

resulting  in  shear  and/or  membrane  locking. 
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For  a  successful  discretization  of  the  theory,  a  consistent  resolution  of  the  locking  deficien¬ 
cies  must  be  sought.  In  Reference  4,  we  have  elaborated  on  an  approximation  strategy  deal¬ 
ing  effectively  with  the  aforementioned  difficulties,  which  involves  a  redefinition  of  the 
penalty  parameters  to  allow  relaxation  of  Equation  14  at  the  element  level  and  an  implementa¬ 
tion  of  appropriate  interpolation  schemes  which  may  best  accommodate  Equation  14.  We 
shall  pursue  both  of  these  avenues  in  deriving  our  three-node  shallow  shell  element. 

Review  ol  Penalty  Relaxation  Concept 

The  first  facet  of  our  approximating  strategy  deals  with  the  issue  of  relaxing  the  enforce¬ 
ment  of  penalty  constraints  at  the  element  level. 

Concurrently  with  the  element  displacement  approximations  denoted  as  w\  uh.  and  9h 
(henceforth,  h  signifies  a  characteristic  length  scale  of  the  discretization),  we  approximate  the 
constitutive  matrices  A  and  G  by  incorporating  the  “element  penalty  relaxation"  parameters: 

N h  «  Aeh,  Qh  ®  **  Gyh,  M*1®  Dich  (15) 

m  s 

where  the  element  strains  are 

eh*  Lx(uh)  -  Lj(O0h*  Yh=  L2(wh)  +  I0h,  ich=  Lj(6h)  (16) 


and  the  penalty  relaxation  parameters  are  nondimensional  positive  quantities  of  the  form: 


(1  +  C^cO  *  (i=  m,s)  (17) 


where  Q  are  positive  element  constants,  and  ax  are  element  analytic  penalty  parameters. 

Note  that  as  =  0(h2/t2)  and  am  =  0[(tf|hh2)2/t2],  where  /c^h  represents  some  characteristic  ini¬ 
tial  curvature  of  the  element.  The  corresponding  principle  of  virtual  work  for  a  single  ele¬ 
ment  approximation  takes  the  form: 

[(N^Se*1  +  (Mh)T«Kh+  (Qh)T6yh  -  q6wh]  dA  -  0 


£ 


(18) 


where  integration  extends  over  the  element  reference  plane  with  Ae  denoting  the  element  ref¬ 
erence  area.  The  resulting  element  strain  energy  appears  in  the  basic  form  of  Equation  13, 
except  that  all  quantities  are  superscribed  with  h  (i.e.,  element  approximations);  however,  the 
element  penalty  parameters  take  a  fundamentally  different  form: 


a.  =  a. /( 1  +  C.o.) 
ii  ii 


(i=  ra,s) . 


(19) 


An  heuristic  argument  in  support  of  Equation  19  may  be  stated  as  follows.  As  t-^0. 
some  degree  of  relaxation  of  constraints.  Equation  14,  may  be  expected  as  now  a  \  approach 
finite  penalty  values  Q1,  which  for  low  order  elements  are  of  0(1). 4  These  thin-regime  pen¬ 
alty  numbers  (C,1)  can  be  seen  to  be  elemcnt/interpolation  dependent.  Thus,  higher  order 
kinematic  interpolations  may  necessitate  larger  C,"1  (i.e.,  smaller  C,),  which  will  ensure  stricter 
element-level  enforcement  of  the  vanishing  penalty  constraints.  Equation  14.  Conversely,  low 
order  kinematic  interpolations  may  require  smaller  C,'1  to  allow  for  greater  element-level  relax 
ation  of  Equation  14.  Clearly,  in  the  limiting  case  when  the  element  kinematic  freedoms  are 
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sufficient  to  capture  the  exact  solution  of  the  problem,  C,'1  can  be  very  large  (or  simply, 

C,  =  0),  in  which  case  Equation  19  takes  on  the  analytic  (unrelaxed)  form. 

From  a  practical  computational  perspective,  Q  can  simply  be  selected  once  and  for  all 
(for  a  particular  element)  from  a  rather  limited  series  of  numerical  tests.  As  the  kinematic 
approximations  improve  with  the  h-refinement  (i.e.,  as  h-*0,  /C|h-*0),  ah ■,  approach  their  ana¬ 
lytic  values  oq,  thus  ensuring  convergence  to  the  “true”  solution  both  in  the  constitutive  and 
kinematic  sense.4,8'11,1*  As  will  be  demonstrated  by  our  numerical  examples,  the  penalty 
parameter  of  Equation  19  is  responsible  for  removing  whatever  spurious  constraining  may 
have  existed  in  the  “unrelaxed”  element  kinematics. 

Anisoparametric  Interpolation  Scheme 

A  complementary  means  for  enhancing  element  behavior  is  by  way  of  appropriate  kine¬ 
matic  interpolations,  termed  anisoparametric,13  which  can  best  accommodate  the  requirements 
of  Equation  14.  The  anisoparametric  strategy  suggests  distinctly  different  degree  polynomials 
for  w,  Q,  and  u  to  reflect  the  differences  in  the  order  of  the  differential  operators  L2  and  1 
in  Equation  14a  and,  likewise,  Lt  and  Lt(£)  in  Equation  14b.  (These  contrast  the  well- 
established  isoparametric  interpolations;  identical  degree  polynomials  for  all  kinematic  vari¬ 
ables.)  The  specific  aim  is  to  design  out  the  unwanted  “spurious”  constraint  equations  that 
may  arise  from  Equation  14. 

To  represent  the  bending  part  of  the  shell  element,  we  adopt  the  three-node 
anisoparametric  plate  element,  '14  in  which  8X  and  dy  are  interpolated  linearly,  while  w  is  rep¬ 
resented  by  a  complete  quadratic  polynomial;  throughout  the  formulation,  area-parametric  coor¬ 
dinates  £  =  (£1;  £2>  £3)  are  used  as  a  basic  for  all  interpolations  (refer  to  Figure  2): 

e£=  N(1)9j  (I=x,y),  vh=  N(2)wh  (20) 

where  Nn^  and  are  the  row  vectors  of  linear  and  quadratic  shape  functions,  respectively, 
and 


(e^)T  =  (6jJ,  (wh)T  =  {w^}  (I=x,y;  j=l,2,3;  k=l,...,6) 

are  the  vectors  of  nodal  dof. 


(21) 


Adopting  the  shell  element  of  constant  curvature  [i.e.,  interpolating  £h(£)  parabolical!}’), 
constraints.  Equation  14b,  necessitate  a  complete  10-term  cubic  polynomial  for  the  u  and  v 
displacements: 


h  „(3)  h  h  „(3)  h 
u  =  N  u  ,  v  =  N  v 

where  is  a  row  vector  of  cubic  shape  functions,  and 

(uh)T  =  (uhk),  (vh)T  =  {vhk>  (k=l , . . . , 10) 


(22) 


(23) 


arc  the  vectors  of  nodal  dof. 
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(a)  Triangle  Area-Parametric  Coordinates 


wn  Nodal  Patterns 


E  q  s  (24) 


ib)  w  Initial  and  Constrained  Nodal  Patterns 


uh  &  vh  Nodal  Patterns 
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(c)  uh  and  vn  Initial  and  Constrained  Nodal  Patterns 


Initial 


Constrained 


Figure  2. 

Evidently,  the  anisoparametric  interpolations  produce  the  same  degree  polynomial  represen 
tations  for  the  left-  and  right-hand  sides  of  the  constraint  equations.  Equation  14;  the  condi¬ 
tion  that  is  paramount  to  improving  element  behavior  in  the  vanishing  thickness  regime. 


Edge  Shear  Constraints 

Although  the  initial  wh  rests  on  six  wh  dof  (i.e.,  three  corner  and  three  midedge  dof).  a 
kinematically  consistent  elimination  of  the  midedge  dof  is  possible  a  priori  to  the  element  still 
ness  derivation.  To  obtain  a  three-node  pattern,  wh  can  be  constrained  by  the  one- 
dimensional  edge  constraints: 

(k) 

*«!  -  M  “h(sv  •£<•>  1  =  °  <k-i-2’3>  <--» 

s  L 

where  s  denotes  a  coordinate  running  along  the  klh  edge  of  the  triangular  element  reference 
plane;  and  0hn(s)  is  the  tangential  edge  rotation  which  is  related  to  0Rx(s)  and  0hv(s)  via  an 
orthogonal  transformation.  From  Equation  24,  there  results  three  decoupled  equations  in 
terms  of  the  midedge  wh  dof,  which  give  rise  to  the  constraints: 

wh=  W  wh+  W  8h  +  W  0h  (25) 

c  o  x  x  y  y 


where  Wq  (q  =  o,  x.  y)  are  3x3  transformation  matrices,  and 


(26) 


/  h\T  .  h  ,  /  h\T  r  h-,  /  .  .  «  *>  \ 

(w  )  =  {w  },  (w  )  =  {w.}  (j=l,2,3). 

c  j+3  j 

Upon  substituting  Equation  25  into  Equation  20,  we  obtain  a  three-node  interpolation  lor 
the  transverse  displacement. 

Edge  Membrane  Constraints 


In  the  manner  analogous  to  the  above  dot'  reduction  for  w  ,  one-dimensional  edge  con¬ 
straints  can  be  devised  to  condense  out  the  intraedge  uh  and  vh  dof.  The  following  con¬ 
straint  equations  provide  four  edge-compatible  relations  for  each  edge: 

(k) 


=  0  (k— 1 ,2, 3 ;  p=l , 2 ) 


£  (s) 

(k) 

r  h  rh  nh  1 

,p 

~P 

u,  -  £  0 
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-  s  s  n 

3sP 

y __(s) 
sn 

=  3sP 

v?  -  sh  eh-  sh  eh 

-  s  n  n  s  s 

where  uh(s)  and  vh(s)  are  cubic  displacement  fields  along  and  normal  to  the  k-th  edge,  respec 
lively,  and 


^  =  5h(0.n 
q  q 


(q=  s,n;  k=l , 2 , 3 > 


V  ° 


128) 


are  the  k-th  edge  slopes. 


By  the  use  of  appropriate  orthogonal  transformations.  Equation  27  is  expressed  in  terms 
of  the  shell  element  variables  of  interest,  namely,  uh,  vh,  0h  and  #hv  dof,  and  algebraically 

L  |.  *  ’  J  ’ 

solved  for  the  intraedge  u  and  v  dof: 

uh  -  u  uh  +  u  eh  +  u  eh 

c  o  xx  y  y 


vh  =  V  vh  +  V  0h  +  V  eh 

c  o  xx  y  y 

where 


,  hxT  ,  h,  ,  hxT  ,  h,  ,  .  ,  nX 

(u)  ={u.},  (v)  ={v.)  ( i=4 , . . . , 9 ) 
c  i  c  1 


(29) 


(30) 


and  Uq  and  Vq  (q  =  o,  x,  y)  are  6x3  transformation  matrices.  Equation  29  is  substituted 
into  the  initial  interpolations.  Equation  22,  to  give  the  constrained  fields  for  the  membrane 
displacements  in  terms  of  the  corner-node  dof  and  two  centroidal  dof.  The  la' ter  dof  are  con 
densed  out  statically  after  the  formation  of  the  element  stiffness  matrix  and  consistent  load 
vector.  Consequently,  a  three-node,  15  dof  element  pattern  is  achieved. 


Note  that  the  edge  constraint  procedures  just  described  preserve  the  original  polynomial 
order  of  the  constrained  variables  (wh,  uh,  and  vh);  moreover,  one  can  show  that  the  con¬ 
strained  fields  are  fully  compatible  across  element  edges,  and  they  allow  for  rigid  body  motion 
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without  straining.  For  further  details  on  this  procedure  and  for  the  explicit  form  of  the 
shape  functions,  refer  to  references  13  and  14. 


The  remainder  of  the  formulation  follows  standard  finite  element  procedures.  Application 
of  the  virtual  work  statement,  Equation  18,  while  performing  exact  integration  (normal  quadra¬ 
ture  rule,1)  throughout,  yields  the  element  stiffness  equations.  The  issue  of  the  rotational  vari¬ 
able  normal  to  the  reference  plane,  #hz,  needed  to  prevent  mathematical  singularities  in  the 
global  coordinates,  produces  three  additional  dof  for  the  element  (e.g.,  see  Reference  10). 

NUMERICAL  EXAMPLES 

An  important  step  in  completing  the  relaxation  methodology  of  the  Finite  Element  Issues 
Section  is  to  obtain  appropriate  a,  parameters  and  the  values  for  C,  (i  =  m,  s).  Herein,  we 
adopt  the  approach  developed  in  Reference  13,  where  as  are  defined  as: 

=  I  k®/  I  k^  (i=s,m;  b  -  bending)  (31) 

in  which  k^  and  k^b  denote  tb element  diagonal  stiffness  coefficients  associated  with  6hx  and 
0hy  dof  for  the  unrelaxed  case,  i.e.,  <p2{  =  1.  As  far  as  the  “optimal”  values  for  Cs  and  Cm, 
these  are  determined  from  numerical  testing.  The  shear  relaxation  constant,  Cs  =  2,  has 
already  been  established  to  ensure  free  of  locking  plate-element  behavior;13  Cm  =  1  was  cho¬ 
sen  from  the  numerical  results  of  the  present  study. 

"M 

The  present  element  was  tested  on  a  series  of  challenging  thin  shell  problems,"  and  a 
moderately  thin  shell.  To  ascertain  the  influence  of  the  membrane  penalty  relaxation  and  the 
membrane  anisoparametric  kinematic  field  upon  element  behavior,  three  versions  of  the  ele¬ 
ment  were  tested.  They  are: 

•  MIN3:  A  facet  triangle  (*£h  =  0)  with  the  shear  relaxation  (Cs  =  2),  possessing  con¬ 
stant  membrane  strains,  no  membrane-bending  coupling,  and  no  membrane  relaxation 
{<P  m  —  1)- 

•  MIN3S:  A  curved  triargle  with  the  shear  relaxation  only  (Cs  =  2,  Cm  =  0). 

•  MIN3SM:  A  curved  element  with  both  the  shear  and  membrane  relaxations  (Cs  =  2. 
Cm  =  1). 

Our  findings  are  summarized  as  follows. 

Test  of  Rigid-Body  Motion:  A  spectral  analysis  was  performed  on  the  element  stiffness 
matrix  for  the  facet,  singly  curved,  and  doubly  curved  element  geometry,  to  check  MIN3’s  abil¬ 
ity  to  move  as  a  rigid  body  without  incurring  any  straining.  Under  all  conditions  tested,  there 
resulted  six  requisite  zero  eigenvalues  associated  with  rigid  body  motion. 

Thin  Cantilevered  Arch:  A  simple  test  of  both  membrane  inextensibility  and  shearless 
deformation  is  a  90°  thin  circular  arch  clamped  at  one  end  and  loaded  by  a  bending  moment 
at  the  other  (Figure  3).  An  additional  modeling  difficulty  here  is  that  the  arch  is  rather  nar¬ 
row  (radius/thickness  ratio:  R/t  =  272;  width/thickness  ratio:  B/t  =  1),  hence  the  element 
aspect  ratios  arc  large.  At  all  discretization  levels,  a  constant  value  of  the  applied  bending 
moment  is  recovered  in  each  element,  with  all  other  stress  resultants  vanishing.  Figure  4 
depicts  a  convergence  study  of  the  tip  bending  rotation,  which  is  also  a  direct  measure  of  the 
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strain  energy  for  this  problem.  Note  that  while  MIN3S  exhibits  considerable  membrane  stiffen¬ 
ing,  MIN3  and  MIN3SM  experience  no  such  difficulty,  converging  rather  rapidly,  with  the 
facet  element  being  slightly  more  flexible.  Both  MIN3  and  MIN3SM  yield  results  of  accept¬ 
able  engineering  accuracy,  even  under  coarse  discretizations. 

CANTILEVERED  ARCH 


NUMBER  OF  NOOES  PER  StOE 


M I N3 
-*-f1IN3S 
-*-f1IN3Sf1 


Figure  4.  90°  thin  cantilevered  arch;  convergence 
study  of  tip  rotation  under  applied  moment. 


Open-Lnded  Pinched  Cylinder:  The  open-ended  cylindrical  shell  subjected  to  two  radial 
forces  ISO  degrees  apart  (Figure  5)  is  a  widely  used  test  problem  to  establish  how  well  a  sin¬ 
gly  curved  shell  element  can  represent  inextensional  bending.1'5  As  t/R-»0,  pure  inextensional 
state  of  deformation  is  attained  in  the  cylinder. 

Figure  6  depicts  a  nine-node-per-side  discretization  (128  elements)  and  its  deformed  shape 
for  a  symmetric  octant  of  the  cylinder.  Figures  7  and  8  show  convergence  studies  of  the 
deflection  under  the  load  for  the  moderately  thin  (R/t  =  50)  and  thin  (R/t  =  2000)  cylinders. 
The  present  results  a'.  . ompared  with  exact  solutions  and  those  obtained  with  three  reduced 
and  one  fully  integ:  Reissner-Mindlin  quadrilateral  elements.  These  elements  are: 

•  9H-S2:  Nine-nod**  eterosis  element  with  selective  integration. 
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•  8S-U2:  Eight-node  serendipity  element  with  uniform  integration. 

•  4L-S1:  Four-node  Lagrange  dement  with  selective  integration. 

•  16L-U4:  Sixteen-node  T  agrangc  fully  integrated  element. 

For  specific  details  on  these  elements  and  references  thereof,  refer  to  References  1  and  5. 


PINCHED  CYLINDER  HTH  OPEN  DCS 


10.35 

4.953  6 
10.5  X  10 
0.3125 


Figure  5.  Pinched,  open-ended  circular  cylinder. 


P/4 


Figure  6.  Pinched,  open-ended  circular  cylinder; 
triangular  element  discretization  for  symmetric 
octant.  Undeformed  and  deformed  shapes. 
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Figure  7.  Moderately  thin-pinched  cylinder  (R/t  =  50);  convergence 
study  of  tip  displacement  under  applied  force. 


Figure  8.  Thin-pinched  cylinder  (R/t  =  20001;  convergence 
study  of  tip  displacement  under  applied  force. 

Clearly,  all  three  MIN3  versions  perform  very  well,  with  MIN3  being  the  stillest  of  the 
three  while  MIN3SM  exhibits  the  most  flexible  behavior.  Evidently,  MIN3SM  is  equally  com¬ 
petitive  with  the  host  performing  quadrilateral  (9H-S2)  for  this  problem. 

Scordelis-Lo  Roof:  The  geometric  description,  loading,  and  material  data  for  this  mem¬ 
brane-response  dominated  cylindrical  shell  are  presented  in  Figure  9.“b  In  Figure  10.  a 
thirteen-node-per-side  mesh  and  the  corresponding  deformed  shape  are  shown  for  a  symmetric- 
quadrant  of  the  shell.  Figure  11  shows  a  convergence  study  of  the  vertical  displacement  at 
the  midpoint  of  the  free  edge  of  the  cylinder.  The  finite  element  displacement  is  normalized 
with  respect  to  a  converged  numerical  solution,  V  =  0.3024. "5 

The  results  obtained  with  four  different  membrane-strain  approximations  for  the  Discrete 
Kirchhoff  Triangle  (DKT)  element"7  (taken  from  Reference  25),  are  also  included  in  this 
study.  Note  that  DKT  is  strictly  a  thin  bending  element  (no  shear  deformation  included). 
Briefly,  the  four  DKT  elements  are: 
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•  DKT/CST:  A  facet  triangle  (no  membrane/bending  coupling;  with  constant  membrane 
strains. 

•  DKT/CST*:  A  constant-strain  triangle  but  which  includes  membrane/bending  coupling 
via  a  membrane  projection  scheme. 

•  DKT/LST:  A  facet  triangle  (no  membrane/bending  coupling)  which  incorporates  a  lin¬ 
ear  membrane  field  and  reduced  integration  of  the  membrane  strain  energy."9 

•  DKT/OB:  A  facet  triangle  (no  membrane/bending  coupling)  with  a  linear  membrane 
field.30 

The  best  performing  elements  for  this  problem  are  DKT/OB,  DKT/LST,  and  MIN3SN1; 
DL'T'/CST  and  DKT/CST*  exhibit  considerable  membrane  stiffening,  showing  a  very  slow  rate 
of  convergence.  Also  note  that  MIN3  and  MIN3S,  although  somewhat  stiller  than  MIN3SM. 
perform  comparably  and  converge  rapidly. 

SCOROELIS-LQ  ROOF 


Figure  9.  Scordeii»-Lo  roof;  vaulted  roof  under  dead  uniform 
load  supported  by  rigid  diaphragms. 


Figure  10,  Scordelis-Lo  roof;  triangular  element  discretization  for 
symmetric  quadrant.  Undeformed  and  deformed  shapes 
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Figure  1 1 .  Scordelis-Lo  roof;  convergence  study  of  vertical 
displacement  at  midpoint  of  free  edge. 


Pinched  Hemisphere:  A  thin  hemispherical  shell  under  self-equilibrating  radial  forces  is  a 
rather  challenging  test  problem  for  doubly  curved  shell  elements.2  The  shell  is  in  the  state 
of  near  extensional  bending,  having  large  rigid-body  rotations  in  the  deformed  configuration. 
Many  commonly  well-behaved  elements,  both  of  quadrilateral  and  triangular  shapes,  exhibit  sig¬ 
nificant  membrane  stiffening  when  modeling  this  problem.5,25 

Figure  12  depicts  a  nine-node-per-side  discretization  of  a  symmetric  quadrant  of  the  hemi¬ 
sphere.  Figure  13  shows  a  convergence  study  for  the  radial  displacement  under  the  applied 
force  which  is  normalized  with  respect  to  an  analytic  solution.31  For  comparison,  results  with 
the  four  DKT  elements  are  also  included.  It  is  evident  that  both  DKT/OB  and  MIN3S 
exhibit  excessive  membrane  stiffening,  while  the  other  elements  suffer  no  such  deficiency  and 
converge  rapidly;  MIN3SM  again  evolves  as  a  reliable  performer. 


PINCHED  HEMISPHERE 

FIXED 


Figure  12.  Pinched  hemisphere;  triangular  element  discretization 
for  symmetric  quadrant. 
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Figure  13.  Pinched  hemisphere;  convergence  study  of 
displacement  under  applied  force. 


Semicircular  Moderately  Thick  Arch:  A  moderately  thick  semicircular  arch  (R/t  =8), 
clamped  at  both  ends  and  subject  to  a  central  vertical  force  (Figure  14),  is  known  to  exhibit 
significant  membrane  and  shear  deformations.11  The  arch  can  be  analyzed,  for  example,  using 
meshes  for  the  pinched  cylinder  problem,  with  the  appropriate  boundary  conditions  and  load¬ 
ing  as  shown  in  Figure  14.  Additionally,  to  ensure  the  x-y  plane  deformations  only,  all  out-of¬ 
plane  global  degrees  of  freedom  are  restrained  in  the  model  (i.e.,  w  =  9X  =  8y  -  0).  The 
analytic  solution  for  the  deflection  under  the  load  can  be  found  in  a  straightforward  manner 
according  to  the  Timoshenko  beam  theory,  and  it  will  serve  as  a  benchmark  in  studying  con¬ 
vergence  for  this  problem.  Figure  15  shows  convergence  results  for  the  deflection  under  the 
load,  which  is  normalized  with  respect  to  the  thin  (inextensional)  solution,  V,hin.  Also  shown 
is  the  exact  solution,  V^ct  =  1.2015  Vthin  (which  includes  shear  and  membrane  deforma¬ 
tions),  and  the  extensional  solution,  Votact  =  1-1299  V,hin  (which  ignores  shear  deformation). 
Note  that  all  three  element  versions  converge  to  the  exact  solution,  with  MIN3SM  yielding 
the  most  accurate  results. 


Figure  14,  Semicircular  moderately  thick  arch  under  vertical 
force  (B/t  =  8)  Undeformed  and  deformed  shapes 
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Figure  15.  Semicircular  moderately  thick  arch;  convergence  study 
of  vertical  displacement  under  applied  force. 

CONCLUSIONS 

A  three-node,  constant  curvature  shallow  shell  element  (MIN3SM)  has  been  formulated 
by  way  of  properly  merging  the  Reissner-Mindlin  theory  of  shear-deformable  plates  and 
Marguerre’s  shallow  shell  equations  for  the  membrane  strains.  The  issues  of  thin-regime 
shear  and  membrane  locking  have  been  addressed  via  shear  and  membrane  penalty  relaxation 
parameters  and  edge-constrained  anisoparametric  interpolations.  The  element  is  fully  inte¬ 
grated  and,  hence,  kinematically  reliable;  it  possesses  six  requisite  rigid-body  modes  and  has 
no  “spurious”  zero-energy  modes. 

Solutions  to  several  locking-sensitive  singly  and  doubly  curved  thin  shell  problems  have 
demonstrated  the  element’s  excellent  modeling  capabilities,  devoid  of  both  shear  and  mem¬ 
brane  locking.  Our  numerical  tests  have  shown  that  MIN3SM  is  consistently  competitive  with 
the  best  DKT  elements.  However,  MIN3SM  has  a  practical  advantage  of  having  a  wider 
range  of  applicability  which  extends  to  moderately  thick  shells. 
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